### Religion Data

library(readxl)
Religion.Dataset <- read_excel("The Religion and State Project, Main Dataset and Societal Module, Round 3.XLSX")

The_Religion_and_State_Project_Main_Dataset_and_Societal_Module_Round_3 -> Religion.Dataset
rm(The_Religion_and_State_Project_Main_Dataset_and_Societal_Module_Round_3)

## Dataset creation year 2014
Main <- dplyr::select(Religion.Dataset, COUNTRY, COWCODE, NXX2014, MXX2014, LXX2014, WSOCDISX2014)
Main2 <- dplyr::select(Religion.Dataset, COUNTRY, COWCODE, NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X)

## Dataset creation year 2000
Main2000 <- dplyr::select(Religion.Dataset, COUNTRY, COWCODE, NXX2000, MXX2000, LXX2000, WSOCDISX2000)


## Building a Dataset by category

###### MX 2014
## MX 01-12
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset, 
                                         MX.01.12 = 
                                           MX01X2014X + 
                                           MX02X2014X + 
                                           MX03X2014X +
                                           MX04X2014X +
                                           MX05X2014X +
                                           MX06X2014X +
                                           MX07X2014X +
                                           MX08X2014X +
                                           MX09X2014X +
                                           MX10X2014X +
                                           MX11X2014X +
                                           MX12X2014X) 

## MX 13-20
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         MX.13.20 = 
                                           MX13X2014X + 
                                           MX14X2014X + 
                                           MX15X2014X +
                                           MX16X2014X +
                                           MX17X2014X +
                                           MX18X2014X +
                                           MX19X2014X +
                                           MX20X2014X) 

## MX 21-27
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         MX.21.27 = 
                                           MX21X2014X + 
                                           MX22X2014X + 
                                           MX23X2014X +
                                           MX24X2014X +
                                           MX25X2014X +
                                           MX26X2014X +
                                           MX27X2014X) 

## MX 28-36
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         MX.28.36 = 
                                           MX28X2014X + 
                                           MX29X2014X + 
                                           MX30X2014X +
                                           MX31X2014X +
                                           MX32X2014X +
                                           MX33X2014X +
                                           MX34X2014X +
                                           MX35X2014X +
                                           MX36X2014X) 

#### NX 2014

## NX 01-05
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         NX.01.05 = 
                                           NX01X2014X + 
                                           NX02X2014X + 
                                           NX03X2014X +
                                           NX04X2014X +
                                           NX05X2014X)


## NX 06-14
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         NX.06.14 = 
                                           NX06X2014X + 
                                           NX07X2014X + 
                                           NX08X2014X +
                                           NX09X2014X +
                                           NX10X2014X +
                                           NX11X2014X +
                                           NX12X2014X +
                                           NX13X2014X +
                                           NX14X2014X)

## NX 15-21
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         NX.15.21 = 
                                           NX15X2014X + 
                                           NX16X2014X + 
                                           NX17X2014X +
                                           NX18X2014X +
                                           NX19X2014X +
                                           NX20X2014X +
                                           NX21X2014X) 

## NX 22-29
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         NX.22.29 = 
                                           NX22X2014X + 
                                           NX23X2014X + 
                                           NX24X2014X +
                                           NX25X2014X +
                                           NX26X2014X +
                                           NX27X2014X +
                                           NX28X2014X +
                                           NX29X2014X) 


#### LX 2014

## LX 01-07                                           
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.01.07 = 
                                           LX01X2014X + 
                                           LX02X2014X + 
                                           LX03X2014X +
                                           LX04X2014X +
                                           LX05X2014X +
                                           LX06X2014X +
                                           LX07X2014X)

## LX 08-11
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.08.11 = 
                                           LX08X2014X + 
                                           LX09X2014X + 
                                           LX10X2014X +
                                           LX11X2014X)

## LX 12-21
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.12.21 = 
                                           LX12X2014X + 
                                           LX13X2014X + 
                                           LX14X2014X +
                                           LX15X2014X +
                                           LX16X2014X +
                                           LX17X2014X +
                                           LX18X2014X + 
                                           LX19X2014X +
                                           LX20X2014X +
                                           LX21X2014X)

## LX 22-26
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.22.26 = 
                                           LX22X2014X + 
                                           LX23X2014X + 
                                           LX24X2014X +
                                           LX25X2014X +
                                           LX26X2014X)

## LX 27-37
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.27.37 = 
                                           LX27X2014X + 
                                           LX28X2014X + 
                                           LX29X2014X +
                                           LX30X2014X +
                                           LX31X2014X +
                                           LX32X2014X + 
                                           LX33X2014X + 
                                           LX34X2014X +
                                           LX35X2014X +
                                           LX36X2014X +
                                           LX37X2014X)

## LX 38-43
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.38.43 = 
                                           LX38X2014X + 
                                           LX39X2014X + 
                                           LX40X2014X +
                                           LX41X2014X +
                                           LX42X2014X +
                                           LX43X2014X)

## LX 44-52
Religion.Dataset.catego <- dplyr::mutate(Religion.Dataset.catego, 
                                         LX.44.52 = 
                                           LX44X2014X + 
                                           LX45X2014X + 
                                           LX46X2014X +
                                           LX47X2014X +
                                           LX48X2014X +
                                           LX49X2014X +
                                           LX50X2014X +
                                           LX51X2014X +
                                           LX52X2014X)




##+#+#+#+##

### 15 VAR Version 2000
###### MX 2000
## MX 01-12
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset, 
                                         MX.01.12 = 
                                           MX01X2000 + 
                                           MX02X2000 + 
                                           MX03X2000 +
                                           MX04X2000 +
                                           MX05X2000 +
                                           MX06X2000 +
                                           MX07X2000 +
                                           MX08X2000 +
                                           MX09X2000 +
                                           MX10X2000 +
                                           MX11X2000 +
                                           MX12X2000) 

## MX 13-20
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         MX.13.20 = 
                                           MX13X2000 + 
                                           MX14X2000 + 
                                           MX15X2000 +
                                           MX16X2000 +
                                           MX17X2000 +
                                           MX18X2000 +
                                           MX19X2000 +
                                           MX20X2000) 

## MX 21-27
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         MX.21.27 = 
                                           MX21X2000 + 
                                           MX22X2000 + 
                                           MX23X2000 +
                                           MX24X2000 +
                                           MX25X2000 +
                                           MX26X2000 +
                                           MX27X2000) 

## MX 28-36
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         MX.28.36 = 
                                           MX28X2000 + 
                                           MX29X2000 + 
                                           MX30X2000 +
                                           MX31X2000 +
                                           MX32X2000 +
                                           MX33X2000 +
                                           MX34X2000 +
                                           MX35X2000 +
                                           MX36X2000) 

#### NX 2000

## NX 01-05
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         NX.01.05 = 
                                           NX01X2000 + 
                                           NX02X2000 + 
                                           NX03X2000 +
                                           NX04X2000 +
                                           NX05X2000)


## NX 06-14
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         NX.06.14 = 
                                           NX06X2000 + 
                                           NX07X2000 + 
                                           NX08X2000 +
                                           NX09X2000 +
                                           NX10X2000 +
                                           NX11X2000 +
                                           NX12X2000 +
                                           NX13X2000 +
                                           NX14X2000)

## NX 15-21
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         NX.15.21 = 
                                           NX15X2000 + 
                                           NX16X2000 + 
                                           NX17X2000 +
                                           NX18X2000 +
                                           NX19X2000 +
                                           NX20X2000 +
                                           NX21X2000) 

## NX 22-29
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         NX.22.29 = 
                                           NX22X2000 + 
                                           NX23X2000 + 
                                           NX24X2000 +
                                           NX25X2000 +
                                           NX26X2000 +
                                           NX27X2000 +
                                           NX28X2000 +
                                           NX29X2000) 


#### LX 2000

## LX 01-07                                           
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.01.07 = 
                                           LX01X2000 + 
                                           LX02X2000 + 
                                           LX03X2000 +
                                           LX04X2000 +
                                           LX05X2000 +
                                           LX06X2000 +
                                           LX07X2000)

## LX 08-11
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.08.11 = 
                                           LX08X2000 + 
                                           LX09X2000 + 
                                           LX10X2000 +
                                           LX11X2000)

## LX 12-21
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.12.21 = 
                                           LX12X2000 + 
                                           LX13X2000 + 
                                           LX14X2000 +
                                           LX15X2000 +
                                           LX16X2000 +
                                           LX17X2000 +
                                           LX18X2000 + 
                                           LX19X2000 +
                                           LX20X2000 +
                                           LX21X2000)

## LX 22-26
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.22.26 = 
                                           LX22X2000 + 
                                           LX23X2000 + 
                                           LX24X2000 +
                                           LX25X2000 +
                                           LX26X2000)

## LX 27-37
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.27.37 = 
                                           LX27X2000 + 
                                           LX28X2000 + 
                                           LX29X2000 +
                                           LX30X2000 +
                                           LX31X2000 +
                                           LX32X2000 + 
                                           LX33X2000 + 
                                           LX34X2000 +
                                           LX35X2000 +
                                           LX36X2000 +
                                           LX37X2000)

## LX 38-43
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.38.43 = 
                                           LX38X2000 + 
                                           LX39X2000 + 
                                           LX40X2000 +
                                           LX41X2000 +
                                           LX42X2000 +
                                           LX43X2000)

## LX 44-52
Religion.Dataset.2000.catego <- dplyr::mutate(Religion.Dataset.2000.catego, 
                                         LX.44.52 = 
                                           LX44X2000 + 
                                           LX45X2000 + 
                                           LX46X2000 +
                                           LX47X2000 +
                                           LX48X2000 +
                                           LX49X2000 +
                                           LX50X2000 +
                                           LX51X2000 +
                                           LX52X2000)


###
Extended1 <- dplyr::select(Religion.Dataset.catego, COUNTRY, COWCODE,
                           MX.01.12:LX.44.52)
Extended2.2000 <- dplyr::select(Religion.Dataset.2000.catego, COUNTRY, COWCODE,
                                MX.01.12:LX.44.52)

Ex2 <- Extended1  ### Backup pathway
Ex2[42, 2] <- NA # Delete Cyprus (Greek)
Ex2[140, 2] <- 345 # Add Serbia
Ex2[158, 2] <- NA # Delete Taiwan
Ex2[182, 2] <- NA # Delete Zanzibar  
  
#######
  ### Main2 Data Editing (Fixing Countries, deleting not official states)
Main3 <- Main2  ### Backup pathway
Main3[42, 2] <- NA # Delete Cyprus (Greek)
Main3[140, 2] <- 345 # Add Serbia
Main3[158, 2] <- NA # Delete Taiwan
Main3[182, 2] <- NA # Delete Zanzibar


#########
## Getting individiual variables
Test <- dplyr::select(Religion.Dataset, COUNTRY, COWCODE, NX01X2014)

## Getting dimension
Test <- dplyr::select(Religion.Dataset, COUNTRY, COWCODE, NXX2014, MXX2014)



# hierachical clustering ####
library(factoextra)
library(fpc)
library(NbClust)

##### TESTS
rel.test1 <- eclust(scale(Main2[,3:6]), "hclust", k = 6, hc_metric = "euclidean", 
                      hc_method = "ward.D2", graph = FALSE, stand = TRUE)

rel.test1$labels <- Main2$COUNTRY

fviz_dend(rel.test1, show_labels = TRUE, type = "rectangle", labels_track_height = 3.35, cex=1.1, lwd=1.4, palette = "jco", main="Title", ylab = "Variables: Spending, Police, Prisonlength, Alternatives", as.ggplot = TRUE)


rownames(rel.test1$data) <- Religion.Dataset$COUNTRY
fviz_cluster(rel.test1, choose.vars=c(1:4) , geom = c("point","text"), 
             ellipse.type = "norm", show.clust.cent = FALSE, repel=TRUE,
             palette = "rainbow", main="Environmental Footprint vs. No. of endangered Species", ggtheme = theme_minimal())
#+#+#+#+#+#+#

### Select Europe
Europe <- dplyr::filter(Main3, COWCODE >= 200 & COWCODE <=399)
### Select Africa
Europe <- dplyr::filter(Main3, COWCODE >= 400 & COWCODE <=626)
### Select Americas
Europe <- dplyr::filter(Main3, COWCODE >= 1 & COWCODE <=180)
### Select Asia
Europe <- dplyr::filter(Main3, COWCODE >= 629 & COWCODE <=818)
### Select SE Asia and Oceania
Europe <- dplyr::filter(Main3, COWCODE >= 819 & COWCODE <=992)



Euro.clust <- eclust(scale(Europe[,3:6]), "hclust", k = 5, hc_metric = "euclidean", 
                     hc_method = "ward.D2", graph = FALSE, stand = TRUE)

Euro.clust$labels <- Europe$COUNTRY

fviz_dend(Euro.clust, show_labels = TRUE, type = "rectangle", labels_track_height = 3.75, cex=1.1, lwd=1.4, palette = "Dark2", main="Religion and State in the Americas", ylab = "Variables: NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X", as.ggplot = TRUE)

## Cluster Stats
Euro.clust$cluster -> Europe$Cluster
knitr::kable(Europe[,c(1,3,4,5,6,7)])
by(Europe, INDICES=list(Europe$Cluster), FUN=summary)


## Mean per Cluster
aggregate(Europe, by=list(Europe$Cluster), FUN=mean, na.rm=TRUE)

## SD per Cluster
aggregate(Europe, by=list(Europe$Cluster), FUN=sd, na.rm=TRUE)


## Analyze 2 dimensional distribution of variables
rownames(Euro.clust$data) <- Europe$COUNTRY


###PCA ####
#+#+#+#+#+###+####+###

library(FactoMineR)
library(corrplot)

PCA.Results <- PCA(Europe[,3:6], scale.unit = TRUE, ncp=5, graph=TRUE)
print(PCA.Results)

## Eigenvalues for each Dimension
get_eigenvalue(PCA.Results)

## Scree plot
fviz_eig(PCA.Results)
fviz_eig(PCA.Results, addlabels = TRUE, ylim = c(0, 50)) # with percentage

##
get_pca_ind(PCA.Results)

### Graph of Variables
PCA.values <- get_pca_var(PCA.Results)
PCA.values$coord

## Nice Corrplot showing the variable importance in Dimensions
corrplot(PCA.values$cos2, is.corr=FALSE)



## Cases (like in fviz_cluster)
fviz_pca_ind(PCA.Results)

## Variables (directions for fviz_cluster)
fviz_pca_var(PCA.Results)
# VAriables NICE WITH COLOR (according to importance!)
fviz_pca_var(PCA.Results, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel= TRUE, arrowsize=1.3, labelsize = 6)




fviz_cluster(Euro.clust, choose.vars=1:4 , geom = c("point","text"), 
             ellipse.type = "norm", show.clust.cent = FALSE, repel=TRUE,
             palette = "rainbow", main="PCA of all 4 Variables included", ggtheme = theme_minimal())





#+#+#+#+#+# not used ####
##
require(devtools)
install_version("backports", version = "1.1.0")

### 1.NbClust
### 2a. Try different regions (Europe) 

## DONE VIA MARKDOWN

### 2b. Try different variable combination for the dimensions

#+#+#+#+#+#+#+#+#+#
### 3a. Democracies & Autocracies (VDEM is included; threshold) #####

`V-Dem-CY-Core-v10` <- readRDS("V-Dem-CY-Core-v10.rds")

library(dplyr)
VDEM.sub <- select(`V-Dem-CY-Core-v10`, country_name, year, COWcode, v2x_polyarchy)
VDEM.sub2 <- filter(VDEM.sub, year == 2014)

VDEM.sub2$dem.aut <- ifelse(VDEM.sub2$v2x_polyarchy >= 0.5, 1, 0)
VDEM.sub2$dem.ano.aut <- ifelse(VDEM.sub2$v2x_polyarchy >= 0.7, 2, 
                                ifelse(VDEM.sub2$v2x_polyarchy >= 0.3, 1, 0))


### Year 2000 (nov 2020 experiment)
VDEM.sub2000 <- filter(VDEM.sub, year == 2000)

VDEM.sub2000$dem.aut <- ifelse(VDEM.sub2000$v2x_polyarchy >= 0.5, 1, 0)
VDEM.sub2000$dem.ano.aut <- ifelse(VDEM.sub2000$v2x_polyarchy >= 0.7, 2, 
                                ifelse(VDEM.sub2000$v2x_polyarchy >= 0.3, 1, 0))

## Colnames for dem.aut
colnames(VDEM.sub2) <- c("name_vdem", "year_vdem", "COWCODE", "vdem_polyarchy", "dem.aut")
colnames(VDEM.sub2000) <- c("name_vdem", "year_vdem", "COWCODE", "vdem_polyarchy", "dem.aut")


## Colnames for dem.ano.aut
colnames(VDEM.sub2) <- c("name_vdem", "year_vdem", "COWCODE", "vdem_polyarchy", "dem.aut", "dem.ano.aut")
colnames(VDEM.sub2000) <- c("name_vdem", "year_vdem", "COWCODE", "vdem_polyarchy", "dem.aut", "dem.ano.aut")


### old code - upgrade with full_join
Main.vdem <- left_join(Main3, VDEM.sub2, by="COWCODE")
Ex2.vdem <- left_join(Ex2, VDEM.sub2, by="COWCODE")

Main2000.vdem <- left_join(Main2000, VDEM.sub2000, by="COWCODE")


#+#+#+#+#+#+#+#+#+#+#+#+#+#+#+#
# DEC 2020 Update
#+#+#


## full_join - dec 2020
library(dplyr)

Main.vdem <- full_join(Main2, VDEM.sub2, by="COWCODE")
Main.vdemB <- full_join(Main2000, VDEM.sub2000, by="COWCODE")

Var15.vdem <- full_join(Extended1, VDEM.sub2, by="COWCODE")
Var15.2000.vdem <- full_join(Extended2.2000, VDEM.sub2000, by="COWCODE")

### 2014 ####

Main.vdem.2014 <- Main.vdem[c(-4,-11,-17,-24,-43:-46,-63:-66,-92,-93:-96,-105,-135:-138,-153:-155,-193:-196,-200,-202,-204),]
Excluded.2014 <- Main.vdem[c(4,11,17,24,43:46,63:66,92,93:96,105,135:138,153:155,193:196,200,202,204),]

Var15.vdem.2014 <- Var15.vdem[c(-4,-11,-17,-24,-43:-46,-63:-66,-92,-93:-96,-105,-135:-138,-153:-155,-193:-196,-200,-202,-204),]


# Recode Serbia because its broken
Main.vdem.2014[130, 7:11]   <- Main.vdem.2014[172, 7:11]
Main.vdem.2014[130, 2]   <- Main.vdem.2014[172, 2]
Main.vdem.2014 <- Main.vdem.2014[-172,]

Main.vdem.2014[169, 1] <- "Dem. Rep. Congo"
Main.vdem.2014 <- Main.vdem.2014[-148,]

## 15 Var
Var15.vdem.2014[130, 18:22]   <- Var15.vdem.2014[172, 18:22]
Var15.vdem.2014[130, 2]   <- Var15.vdem.2014[172, 2]
Var15.vdem.2014 <- Var15.vdem.2014[-172,]

Var15.vdem.2014[169, 1] <- "Dem. Rep. Congo"
Var15.vdem.2014 <- Var15.vdem.2014[-148,]


## writing to file


write.csv(Main.vdem.2014, "Main_vdem_2014.csv", row.names = FALSE) 
write.csv(Excluded.2014, "excluded_2014.csv", row.names = FALSE) 

write.csv(Var15.vdem.2014, "Var15_vdem_2014.csv", row.names = FALSE) 



library(arsenal)
comparedf(Main.vdem.2014, Main.vdem)
summary(comparedf(Main.vdem.2014, Main.vdem))

### 2000 ####
Main.vdem.2000 <- Main.vdemB[c(-4,-11,-17,-24,-43:-45,-62:-64,-90:-93,-102,-117,-132:-134,-149,-150,-159,-168,-172,-188:-190,-194),]
Excluded.2000 <- Main.vdemB[c(4,11,17,24,43:45,62:64,90:93,102,117,132:134,149,150,159,168,172,188:190,194),]

Var15.vdem.2000 <- Var15.2000.vdem[c(-4,-11,-17,-24,-43:-45,-62:-64,-90:-93,-102,-117,-132:-134,-149,-150,-159,-168,-172,-188:-190,-194),]


Main.vdem.2000[129, 7:11]   <- Main.vdem.2000[169, 7:11]
Main.vdem.2000[129, 2]   <- Main.vdem.2000[169, 2]
Main.vdem.2000 <- Main.vdem.2000[-168:-170,]

Main.vdem.2000[165, 1] <- "Dem. Rep. Congo"


Var15.vdem.2000[129, 18:22]   <- Var15.vdem.2000[169, 18:22]
Var15.vdem.2000[129, 2]   <- Var15.vdem.2000[169, 2]
Var15.vdem.2000 <- Var15.vdem.2000[-168:-170,]

Var15.vdem.2000[165, 1] <- "Dem. Rep. Congo"


write.csv(Main.vdem.2000, "Main_vdem_2000.csv", row.names = FALSE) 
write.csv(Excluded.2000, "excluded_2000.csv", row.names = FALSE)

write.csv(Var15.vdem.2000, "Var15_vdem_2000.csv", row.names = FALSE) 

#+#+#+#+#+#+#+#+#+#+#+#+#

na.omit(Main.vdem) -> Main.vdem2
na.omit(Ex2.vdem) -> Ex2.vdem2

na.omit(Main2000.vdem) -> Main2000.vdem2

# year 2014 version

Main.vdem.demo <- filter(Main.vdem2, dem.aut==1)
Main.vdem.auto <- filter(Main.vdem2, dem.aut==0)

Ex2.vdem.demo <- filter(Ex2.vdem2, dem.aut==1)
Ex2.vdem.auto <- filter(Ex2.vdem2, dem.aut==0)

Main.vdem.demo2 <- filter(Main.vdem2, dem.ano.aut==2)
Main.vdem.ano2 <- filter(Main.vdem2, dem.ano.aut==1)
Main.vdem.auto2 <- filter(Main.vdem2, dem.ano.aut==0)

# year 2000 version
Main2000.vdem.demo <- filter(Main2000.vdem2, dem.aut==1)
Main2000.vdem.auto <- filter(Main2000.vdem2, dem.aut==0)

Main2000.vdem.demo2 <- filter(Main2000.vdem2, dem.ano.aut==2)
Main2000.vdem.ano2 <- filter(Main2000.vdem2, dem.ano.aut==1)
Main2000.vdem.auto2 <- filter(Main2000.vdem2, dem.ano.aut==0)


## write for markdown
setwd("C:/Users/Felix Ettensperger/Desktop/ECPR Draft Religion/VDEM Data")

## Main Data
write.csv(Main.vdem.demo2, file="Demo07.csv")
write.csv(Main.vdem.demo, file="Demo05.csv")
write.csv(Main.vdem.auto2, file="Auto03.csv")
write.csv(Main.vdem.auto, file="Auto05.csv")
write.csv(Main.vdem.ano2, file="Ano0307.csv")

## Data with 15 Categories 
write.csv(Ex2.vdem.demo, file="Cat15_Demo05.csv")
write.csv(Ex2.vdem.auto, file="Cat15_Auto05.csv")
##+#+#+#+#+# Version with 15 Categories and VDEM Data

## Data with year 2000 in RaS and VDEM Data
write.csv(Main2000.vdem.demo2, file="y2000_Demo07.csv")
write.csv(Main2000.vdem.demo, file="y2000_Demo05.csv")
write.csv(Main2000.vdem.auto2, file="y2000_Auto03.csv")
write.csv(Main2000.vdem.auto, file="y2000_Auto05.csv")
write.csv(Main2000.vdem.ano2, file="y2000_Ano0307.csv")



###########
## Pew Data Clustering ## Add Demo-Auto
library(readxl)
pewdata <- read_excel("Pew Religion Data/Pew with COW.xlsx")
Dataset <- pewdata
Dataset2 <- pewdata

## Create new variables ## Elinas Script
# Creating the Social Hostilities Index
Dataset$SocialHostilities <- Dataset$SHI_Q_1_HARASSMENT+Dataset$SHI_Q_2+Dataset$SHI_Q_3+Dataset$SHI_Q_4+Dataset$SHI_Q_5+Dataset$SHI_Q_6+Dataset$SHI_Q_7+Dataset$SHI_Q_8+Dataset$SHI_Q_9+Dataset$SHI_Q_10+Dataset$SHI_Q_11+Dataset$SHI_Q_12+Dataset$SHI_Q_13
Dataset2$SocialHostilities <- Dataset$SHI_Q_1_HARASSMENT+Dataset$SHI_Q_2+Dataset$SHI_Q_3+Dataset$SHI_Q_4+Dataset$SHI_Q_5+Dataset$SHI_Q_6+Dataset$SHI_Q_7+Dataset$SHI_Q_8+Dataset$SHI_Q_9+Dataset$SHI_Q_10+Dataset$SHI_Q_11+Dataset$SHI_Q_12+Dataset$SHI_Q_13



# Creating the Government Restrictions Index
Dataset$GovernmentRegulation <- Dataset$GRI_Q_1+Dataset$GRI_Q_2+Dataset$GRI_Q_3+Dataset$GRI_Q_4+Dataset$GRI_Q_5+Dataset$GRI_Q_6+Dataset$GRI_Q_7+Dataset$GRI_Q_8+Dataset$GRI_Q_9+Dataset$GRI_Q_10+Dataset$GRI_Q_11+Dataset$GRI_Q_12+Dataset$GRI_Q_13+Dataset$GRI_Q_14+Dataset$GRI_Q_15+Dataset$GRI_Q_16+Dataset$GRI_Q_17+Dataset$GRI_Q_18+Dataset$GRI_Q_19


# Subindicator in Dataset 2: Government laws and policies restricting religious freedom
Dataset2$Freedom<-Dataset$GRI_Q_1+Dataset$GRI_Q_2+Dataset$GRI_Q_3+Dataset$GRI_Q_14+Dataset$GRI_Q_18
# Subindicator in Dataset 2: Government limits on activities of religious groups and individuals
Dataset2$Groups<-Dataset$GRI_Q_4+Dataset$GRI_Q_5+Dataset$GRI_Q_6+Dataset$GRI_Q_7+Dataset$GRI_Q_8+Dataset$GRI_Q_9+Dataset$GRI_Q_10
# Subindicator in Dataset 2: Government harassment of religious groups
Dataset2$Harassment<-Dataset$GRI_Q_11+Dataset$GRI_Q_12+Dataset$GRI_Q_13+Dataset$GRI_Q_15+Dataset$GRI_Q_16+Dataset$GRI_Q_17+Dataset$GRI_Q_19


# Creating the Government Favoritism Index
Dataset$Favoritism <- Dataset$GRI_Q_20_1+Dataset$GRI_Q_20_2+Dataset$GRI_Q_20_3+Dataset$GRI_Q_20_4+Dataset$GRI_Q_20_5
Dataset2$Favoritism <- Dataset$GRI_Q_20_1+Dataset$GRI_Q_20_2+Dataset$GRI_Q_20_3+Dataset$GRI_Q_20_4+Dataset$GRI_Q_20_5


colnames(Dataset)[5] <- "COWCODE"
colnames(Dataset)[4] <- "COUNTRY"

colnames(Dataset2)[5] <- "COWCODE"
colnames(Dataset2)[4] <- "COUNTRY"

VDEM.sub2$COWCODE <- as.character(VDEM.sub2$COWCODE)
pewdata.vdem <- left_join(Dataset, VDEM.sub2, by="COWCODE")
pewdata.vdem2 <- left_join(Dataset2, VDEM.sub2, by="COWCODE")

pewdata.vdem[,c(2,4,5,84:90)] -> pew.vdem
na.omit(pew.vdem) -> pew.vdem2

pewdata.vdem2[,c(2,4,5,84:92)] -> pew.subcluster
na.omit(pew.subcluster) -> pew.subcluster2



pew.vdem.demo <- filter(pew.vdem2, dem.aut==1)
pew.vdem.auto <- filter(pew.vdem2, dem.aut==0)

pew.sub.demo <- filter(pew.subcluster2, dem.aut==1)
pew.sub.auto <- filter(pew.subcluster2, dem.aut==0)

## write for markdown
setwd("C:/Users/Felix Ettensperger/Desktop/ECPR Draft Religion/VDEM Data")

write.csv(pew.vdem.demo, file="pew_demo05.csv")
write.csv(pew.vdem.auto, file="pew_auto05.csv")

write.csv(pew.sub.demo, file="pew_sub_demo05.csv")
write.csv(pew.sub.auto, file="pew_sub_auto05.csv")






######

## Cluster Dendrograms

## demo >0.5
Dem.clust <- eclust(scale(Main.vdem.demo[,3:6]), "hclust", k = 5, hc_metric = "euclidean", 
                     hc_method = "ward.D2", graph = FALSE, stand = TRUE)

Dem.clust$labels <- Main.vdem.demo$COUNTRY

fviz_dend(Dem.clust, show_labels = TRUE, type = "rectangle", 
          labels_track_height = 3.75, cex=1.1, lwd=1.4, palette = "Dark2", 
          main="Religion and State in Democracies (vdem_poly >0.5)", 
          ylab = "Variables: NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X", 
          as.ggplot = TRUE)

## demo >0.7
Dem.clust2 <- eclust(scale(Main.vdem.demo2[,3:6]), "hclust", k = 5, hc_metric = "euclidean", 
                    hc_method = "ward.D2", graph = FALSE, stand = TRUE)

Dem.clust2$labels <- Main.vdem.demo2$COUNTRY

fviz_dend(Dem.clust2, show_labels = TRUE, type = "rectangle", 
          labels_track_height = 3.75, cex=1.1, lwd=1.4, palette = "Dark2", 
          main="Religion and State in Democracies (vdem_poly >0.7)", 
          ylab = "Variables: NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X", 
          as.ggplot = TRUE)

## auto >0.7
Auto.clust2 <- eclust(scale(Main.vdem.auto2[,3:6]), "hclust", k = 5, hc_metric = "euclidean", 
                     hc_method = "ward.D2", graph = FALSE, stand = TRUE)

Auto.clust2$labels <- Main.vdem.auto2$COUNTRY

fviz_dend(Auto.clust2, show_labels = TRUE, type = "rectangle", 
          labels_track_height = 3.75, cex=1.1, lwd=1.4, palette = "Dark2", 
          main="Religion and State in Autocracies (vdem_poly <0.3)", 
          ylab = "Variables: NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X", 
          as.ggplot = TRUE)


## auto >0.5
Auto.clust <- eclust(scale(Main.vdem.auto[,3:6]), "hclust", k = 5, hc_metric = "euclidean", 
                      hc_method = "ward.D2", graph = FALSE, stand = TRUE)

Auto.clust$labels <- Main.vdem.auto$COUNTRY

fviz_dend(Auto.clust, show_labels = TRUE, type = "rectangle", 
          labels_track_height = 4.75, cex=1.1, lwd=1.4, palette = "Dark2", 
          main="Religion and State in Autocracies (vdem_poly <0.5)", 
          ylab = "Variables: NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X", 
          as.ggplot = TRUE)

## ano >0.3 <0.7
Ano.clust <- eclust(scale(Main.vdem.ano2[,3:6]), "hclust", k = 5, hc_metric = "euclidean", 
                     hc_method = "ward.D2", graph = FALSE, stand = TRUE)

Ano.clust$labels <- Main.vdem.ano2$COUNTRY

fviz_dend(Ano.clust, show_labels = TRUE, type = "rectangle", 
          labels_track_height = 4.75, cex=1.1, lwd=1.4, palette = "Dark2", 
          main="Religion and State in Anocracies (vdem_poly >0.3 & <0.7)", 
          ylab = "Variables: NXX2014X, MXX2014X, LXX2014X, WSOCDISX2014X", 
          as.ggplot = TRUE)


### 3b. Stateness / Rule of Law (Indicator) 
### 2b note: Elina: Different sub-variables impact on overall dimensions

### 3c. Run it with 3 indicators not 4 (WSOCDISX2014)
### 3d. Deconstruct NX into subcategories = 4 Variables (adding nx1-7 to bla etc)
### 3e. Deconstruct LX into subcategories = 7 Variables (adding lx1-7 to bla etc)
### 3f. Use Decostructed Variables LX and NX together (11 Variables) 

##++## Analysis and Results ####
### 5. Average for cluster
### 6. Automatisation of Web-Diagram (Radar Chart)
### 7. consider years
### 8. country change from 1990-2014 (Sankey Diagram)



